#  load genome-wide clones file
load_raw_gw_clones_data <- function ( file ) {
  
  raw_gw_clones_data<- read.table(
    file,
    skip=1,
    colClasses = c(
      "NULL", "character", "character", "NULL", "NULL",
      "NULL", "character", "numeric", "numeric", "NULL"
    ),
    col.names = c(
      "NULL", "word_2nd", "word_3rd", "NULL", "NULL",
      "NULL", "contig_id", "start", "end", "NULL"
    ),
    stringsAsFactors =  FALSE 
    
  )
  
  #  remove entries for which word_2nd starts with "NC"
  
  raw_gw_clones_data <- raw_gw_clones_data[ !( substr( raw_gw_clones_data$word_2nd, 1, 2 ) %in% "NC" ), ]
  row.names( raw_gw_clones_data ) <- NULL
  
  
  #  correcting for size2 problem in word_3rd
  
  size2_indices <-  nchar( raw_gw_clones_data$word_3rd ) == 2
  size2_correction <- paste( 
    substr( raw_gw_clones_data$word_3rd[ size2_indices ], 1, 1 ),
    "0",
    substr( raw_gw_clones_data$word_3rd[ size2_indices ], 2, 2 ),
    sep = ""
  )
  raw_gw_clones_data$word_3rd[ size2_indices ] <- size2_correction;
  
  #test cat( "size2 correction done: ", !( FALSE %in% (nchar(raw_gw_clones_data$word_3rd) == 3) ), "\n" )
  
  #  correcting for X and G problem in word_2nd
  
  #  X problem
  char_1st <- substr( raw_gw_clones_data$word_2nd, 1, 1 )
  digits <- as.numeric( substring( raw_gw_clones_data$word_2nd, 2) )
  indices <- char_1st == "X" && digits >= 1 && digits <= 25
  digits[ indices ] <- digits[ indices ] + 106
  
  #  G problem
  indices <- char_1st == "G" && digits >= 1 && digits <= 25
  digits[ indices ] <- digits[ indices ] + 131
  char_1st[ indices ] <- "X"
  
  raw_gw_clones_data$word_2nd <- paste( char_1st, as.character( digits ), sep = "" )
  
  
  #  join word_2nd and word_3rd
  raw_gw_clones_data <- cbind( raw_gw_clones_data, gw_clone_names = paste( raw_gw_clones_data$word_2nd, raw_gw_clones_data$word_3rd, sep = "") )
  raw_gw_clones_data$word_2nd <- NULL
  raw_gw_clones_data$word_3rd <- NULL
  
  #  get contig id
  raw_gw_clones_data$contig_id <- as.numeric( matrix( unlist( strsplit( raw_gw_clones_data$contig_id, "\\." ) ), byrow = TRUE, ncol = 2 )[,2] )
  
  # put the variables in order
  raw_gw_clones_data <- raw_gw_clones_data[c("gw_clone_names","contig_id", "start", "end")]
  
  rm( list=c( "indices", "digits", "char_1st", "size2_indices", "size2_correction" ) )
  
  return (raw_gw_clones_data)
}

#  load genome-wide ncu file
load_raw_gw_ncu_data <- function ( file ) {
  
  raw_gw_ncu_data <- read.table(
    file,
    colClasses = c(
      "NULL", "NULL", "NULL", "character", "NULL",
      "NULL", "numeric", "numeric", "NULL", "NULL",
      "NULL", "NULL", "character", "NULL", "NULL",
      "NULL","NULL"
    ),
    col.names = c(
      "NULL", "NULL", "NULL", "contig_id", "NULL",
      "NULL", "start", "end", "NULL", "NULL",
      "NULL", "NULL", "ncu", "NULL", "NULL",
      "NULL","NULL"
    ),
    quote = "\"",
    stringsAsFactors =  FALSE 
    
  );
  
  #  get contig id
  raw_gw_ncu_data$contig_id <-  matrix( unlist( strsplit( raw_gw_ncu_data$contig_id, "\\." ) ), byrow = TRUE, ncol = 2 )[,2] 
  
  size <- length ( unique(raw_gw_ncu_data$ncu) )
  temp_data <- data.frame( 
    ncu = rep("",size),
    contig_id = rep("",size), 
    start =  rep(0,size), 
    end = rep(0,size), 
    stringsAsFactors =  FALSE 
  )
  
  ncu <- raw_gw_ncu_data[ 1, ]$ncu
  start_pos <- raw_gw_ncu_data[ 1, ]$start
  end_pos <- raw_gw_ncu_data[ 1, ]$end
  contig_id <- raw_gw_ncu_data[ 1, ]$contig_id
  
  count <- 1
  for( i in 1 : nrow( raw_gw_ncu_data ) ) {
    
    record <- raw_gw_ncu_data[ i, ];
    
    if ( record$ncu == ncu ) {
      
      start_pos <- c ( start_pos, record$start )
      
      end_pos <- c ( end_pos, record$end )
      
    } else {
      
      temp_data[ count, ] <- c( contig_id, ncu, min( start_pos ), max ( end_pos ) ) 
      
      start_pos <- raw_gw_ncu_data[ i+1, ]$start
      
      end_pos <- raw_gw_ncu_data[ i+1, ]$end
      
      count <- count + 1;
      
    }
    
    ncu = record$ncu
    contig_id = record$contig_id
    
  }
  
  temp_data <- temp_data[ -1, ]
  row.names( temp_data ) <- NULL	
  
  raw_gw_ncu_data <- temp_data
  raw_gw_ncu_data$start <- as.numeric(raw_gw_ncu_data$start)
  raw_gw_ncu_data$end <- as.numeric(raw_gw_ncu_data$end)
  
  rm( list= c( "temp_data" , "record", "start_pos", "end_pos", "contig_id", "count" ) )
  
  return (raw_gw_ncu_data)
}

#  load gw_gene_ncu.txt file
load_raw_gw_gene_ncu_data <- function ( file ) {
  
  raw_gw_gene_ncu_data <- read.csv(
    header = FALSE,
    file = file,
    colClasses = c(
      "character",
      "character"
    ),
    col.names = c(
      "gene_name",
      "ncu"
    ),
    stringsAsFactors =  FALSE 
  );
  
  return (raw_gw_gene_ncu_data)
}


#  load gw_compl.txt file
load_raw_compl_data <- function ( file ) {
  
  raw_compl_data <- read.table(
    file,		
    fill = FALSE,
    stringsAsFactors =  FALSE ,
  );
  
  temp1 <- vector()
  mylist <- strsplit( raw_compl_data$V1, "\\," )
  for ( i in 1: length( mylist ) ) {
    temp1 <- c(temp1,mylist[[i]][1])
    mylist[[i]] <- mylist[[i]][-1]
  }
  names(mylist) <- temp1
  raw_compl_data <- mylist
  rm( list=c( "temp1", "mylist","i" ) )
  
  return (raw_compl_data)
}

# load clone_names_<i>
load_ch_clone_names_data <- function( file ) {
  
  scan( file, what = character(0) )
  
}

# load probe_names_<i>
load_ch_probe_names_data <- function( file ) {
  
  scan( file, what = character(0) )
  
}

# load hybridization_data_<i>
load_ch_hybridization_data <- function ( file ) {
  read.csv (
    file,
    header = FALSE
  )
}

# load gene_names_<i>
load_ch_gene_names_data <- function( file ) {
  
  scan( file, what = character(0) )
  
}

# load genotype_data_<i>
load_ch_genotype_data <- function ( file ) {
  read.csv (
    file,
    header = FALSE,
    na.strings = "9"
  )
}

check_IC_1 <- function () {
  
  dup1 <- duplicated(raw_gw_gene_ncu_data$gene)
  dup2 <- duplicated(raw_gw_gene_ncu_data$ncu)
  if( TRUE %in% dup1 ){
    cat ("There are duplicate genes in gw_gene_ncu. They are \n",
         paste( raw_gw_gene_ncu_data$gene[ dup1 ], "\n", sep=""), "\n"
    )
    stop( "correct the duplicate genes " )
  }
  if( TRUE %in% dup2 ){
    cat ("There are duplicate ncu in gw_gene_ncu. They are \n",
         paste( raw_gw_gene_ncu_data$ncu[ dup2 ], "\n", sep=""), "\n"
    )
    stop( "correct the duplicate ncu" )
  }
  
  rm( list= c ("dup1","dup2") )
  
}


check_IC_2 <- function () {
  temp1 <- raw_gw_gene_ncu_data$ncu %in% raw_gw_ncu_data$ncu
  if( length(temp1) != length(raw_gw_gene_ncu_data$ncu) ){
    cat ( "Some ncu could not be found in gw_ncu. They are: \n",
          paste ( raw_gw_gene_ncu_data$ncu[ !temp1], "\n", sep = "" ), "\n",
          "The gene-ncu pairs corresponding to these ncu will be removed from gw_gene_ncu \n"
    )
    raw_gw_gene_ncu_data <<- raw_gw_gene_ncu_data[ temp1, ]
  }
  rm( list= c ("temp1") )
}	


check_IC_3 <-  function() {
  temp1 <- raw_gw_ncu_data$ncu %in% raw_gw_gene_ncu_data$ncu
  if( FALSE %in% temp1 ){
    raw_gw_ncu_data <<- raw_gw_ncu_data[ raw_gw_ncu_data$ncu %in% raw_gw_gene_ncu_data$ncu, ]
    row.names( raw_gw_ncu_data ) <<- NULL
  }
  rm( list= c ("temp1") )
}

# Note that raw_gw_compl_data will not be written to disk. It only checks for error on the original source.
# Edit the original for errors and include it manually in the managed-input directory. The reason for doing it automatically is
# 1. It does not need much editing 2. I could not managed to print the list in the format desired.
check_IC_4 <- function () {
  dup1 <- duplicated(names(raw_gw_compl_data))
  if( TRUE %in% dup1 ){
    cat ("There are duplicate genes in gw_compl_data. They are \n",
         paste( names(raw_gw_compl_data)[ dup1 ], "\n", sep=""), "\n"
    )
    stop( "remove the duplicate genes" )
    #raw_gw_compl_data <<- raw_gw_compl_data[ !dup1 ]  #uncomment if you ever want to write to disk automatically.
  }
  for( i in 1: length(raw_gw_compl_data) ){
    dup2 <- duplicated( raw_gw_compl_data[[i]] )
    if ( TRUE %in% dup2 ){
      cat ("There are duplicate cosmids in gw_compl_data for gene ", 
           names(raw_gw_compl_data)[i], 
           " They are \n",
           paste( raw_gw_compl_data[[i]][ dup2 ], "\n", sep=""), "\n"
      )
      rm( list= c ("dup1","i","dup2") )
      stop( "correct the duplicate cosmids!" )
    }
    rm( list= c ("dup2") )
  }
  rm( list= c ("dup1","i") )
}


check_IC_5 <-  function() {
  
  all_clones <- vector()
  
  for(linkage_group in 1:7 ){
    
    # check duplicates from clones
    if( TRUE %in% duplicated(ch_clone_names_data[[ linkage_group ]]) ){
      stop( paste("duplicates in clone names. linkage group: ", linkage_group ))
    }
    
    all_clones <- c( all_clones, ch_clone_names_data[[ linkage_group ]])
    temp1 <- ch_probe_names_data[[ linkage_group ]] %in% ch_clone_names_data[[ linkage_group ]]
    if ( FALSE %in% temp1 )	{
      temp2 <- ch_probe_names_data[[ linkage_group ]][ !temp1];
      cat ( "Some probes for linkage group ", linkage_group, " are not present in ch_clones.\n",
            paste( c("Previous Count: ", " To be removed: "), c( length(temp1), length(temp2) ), sep=""), "They are: \n",
            paste(temp2,"\n",sep=""),
            "These probes will now be removed \n"
      )
      
      ch_probe_names_data[[ linkage_group ]] <<- ch_probe_names_data[[ linkage_group ]][ temp1 ]
      
      # check duplicates from probes
      if( TRUE %in% duplicated(ch_probe_names_data[[ linkage_group ]]) ){
        stop( paste("duplicates in probe names. linkage group: ", linkage_group ))
      }
      
      update_hybridization_data(linkage_group)
      rm( list= c ("temp2") )
    }
  }
  temp2 <- raw_gw_clones_data$gw_clone_names %in% all_clones
  raw_gw_clones_data <<- raw_gw_clones_data[ temp2, ]
  
  # check duplicates from raw_gw_clones_data
  if( TRUE %in% duplicated(raw_gw_clones_data$gw_clone_names) ){
    stop("duplicates in gw clone names.")
  }
  
  rownames(raw_gw_clones_data) <- NULL
  rm( list= c ("temp1","temp2","linkage_group") )
}

check_IC_6 <- function() {
  
  for(linkage_group in 1:7 ){
    
    temp2 <- as.matrix(ch_hybridization_data[[ linkage_group ]][ch_probe_names_data[[ linkage_group ]],])
    
    # find probes that did not hybridize to themselves
    temp3 <- vector()
    for( i in 1:nrow(temp2) ){
      if ( temp2[i,i] != 1) temp3 <- c(temp3, rownames(temp2)[i])
    }
    if( length(temp3) > 0 ){
      cat( length(temp3), "out of", ncol(temp2), "probes did not hybridize to themselves. They are: \n",
           paste(temp3,"\n",sep=""),
           "These probes will now be removed.\n"
      )
    }
    
    black_list <- temp3
    
    # find the overlapping probes
    temp3 <- vector()
    temp4 <- list()
    for( i in 1:nrow(temp2) ){
      j = i+1
      temp5 <- vector()
      while( j <= ncol(temp2) ){
        if ( temp2[i,j] != 0) {
          temp5 <- c(temp5, rownames(temp2)[j])
        }
        j = j+1
      }
      if(length(temp5) > 0 ) {
        temp4[ rownames(temp2)[i] ] <- list(temp5)
        temp3 <- c(temp3, temp5)
      }
      
      rm(list=c("j","temp5"))
    }
    if( length(temp3) > 0 ){
      
      temp3 <- unique(temp3)
      cat( length(temp3),"out of", ncol(temp2), "probes were found to be overlapping. They are: \n") 
      print(temp4)
      
      # remove overlapping probes while trying to keep probes that have max number of gene tags
      for( i in 1:length(temp4) ){
        list1 <- c( names(temp4)[i], temp4[[i]] )
        list1 <- list1[ !(list1 %in% black_list) ]
        list1<- list1[-order(find_tag_count(list1)[1])]
        black_list <- c( black_list, list1)
      }
      
      keep <- !(ch_probe_names_data[[ linkage_group ]] %in% black_list)
      ch_probe_names_data[[ linkage_group ]] <<- ch_probe_names_data[[ linkage_group ]][ keep ]
      update_hybridization_data(linkage_group)
      
      rm(list=c("i","keep","list1"))
    }
  }
  rm(list=c("temp2","temp3","temp4","black_list"))
  
}

check_IC_7 <-  function() {
  for(linkage_group in 1:7 ){
    # keep odd number of genes
    if( length(ch_gene_names_data[[ linkage_group ]]) %% 2 == 0 ){
      ch_gene_names_data[[ linkage_group ]]<<- ch_gene_names_data[[ linkage_group ]][ -1 ]
      update_genotype_data(linkage_group)
    }
  }
}


update_genotype_data <- function(linkage_group) {
  
  temp1 <- colnames(ch_genotype_data[[ linkage_group ]]) %in% ch_gene_names_data[[ linkage_group ]]
  
  ch_genotype_data[[ linkage_group ]] <<- ch_genotype_data[[ linkage_group ]][ temp1,]
  
  rm( list= c ("temp1") )
}

find_tag_count <- function (probes){
  tags <- vector()
  for(i in 1:length(probes)){
    temp1 <- probe_tag_info$probes %in% probes[i]
    if( TRUE %in% temp1 ){
      tags[i] <- probe_tag_info$seq_tag[temp1] + probe_tag_info$compl_tag[temp1]
    } else{
      tags[i] <- 0
    }
  }
  return (tags)
}

update_hybridization_data <- function ( linkage_group ){
  
  temp1 <- colnames(ch_hybridization_data[[ linkage_group ]]) %in% ch_probe_names_data[[ linkage_group ]]
  
  ch_hybridization_data[[ linkage_group ]] <<- ch_hybridization_data[[ linkage_group ]][ ,temp1]
  
  rm( list= c ("temp1") )
}

create_probe_tag_info <- function (){
  probes <- vector()
  for(i in 1:7){
    probes <- c( probes, ch_probe_names_data[[ i ]])
  }
  tagged_probes <- vector()
  for(i in 1:length(raw_gw_compl_data)){
    tagged_probes <- c( tagged_probes, raw_gw_compl_data[[i]])
  }
  probes <- unique(probes)
  probe_tag_info <- data.frame( probes = probes, seq_tag = rep(0,length(probes)), compl_tag = rep(0,length(probes)) )
  
  temp1 <- raw_gw_clones_data[raw_gw_clones_data$gw_clone_names %in% probes, ]
  for( i in 1:nrow(temp1) ){
    temp2 <- raw_gw_ncu_data$contig_id == temp1$contig_id[i]
    temp <- (raw_gw_ncu_data$start[temp2] >= temp1$start[i] & raw_gw_ncu_data$start[temp2] <= temp1$end[i]) |
      (raw_gw_ncu_data$end[temp2] >= temp1$start[i] & raw_gw_ncu_data$end[temp2] <= temp1$end[i])
    
    index <- probe_tag_info$probes %in% temp1$gw_clone_names[i]
    
    probe_tag_info$seq_tag[ index ] = sum(temp)
    probe_tag_info$compl_tag[ index ] = length( grep( temp1$gw_clone_names[i], tagged_probes ) )
  }
  
  for( tagged_probe in tagged_probes){
    probe_tag_info$compl_tag[ probe_tag_info$probes == tagged_probe] = length( grep( tagged_probe, tagged_probes ) )
  }
  
  probe_tag_info <- probe_tag_info[ probe_tag_info$seq_tag > 0 | probe_tag_info$compl_tag > 0,] 
  rownames(probe_tag_info) <- NULL
  
  rm( list= c ("probes", "tagged_probes", "temp1", "temp2", "temp") )
  return( probe_tag_info )
  
}




